library(ChromCom)
library(ggplot2)
library(gridExtra)
library(reshape2)

11/07/2017

Writing introduction. First script to do simulation.

12/07/2017

Let’s try if it works. An example for \(t_1 = -60\) min, \(\Delta t = 0\), \(k_1 = 0.04\ {\rm min}^{-1}\) and \(k_2 = 0.04\ {\rm min}^{-1}\).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0.04
)
chr <- ChromCom3(pars)
chr <- generateCells(chr, nsim=1000)
plotTimelines(chr)

Just to check if output is as expected. With \(k_2=0\) we should have pure exponential decay. The yellow curve is \(e^{-k_1 t}\).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
)
tchr <- ChromCom3(pars)
tchr <- generateCells(tchr, nsim=1000)
x <- seq(from=pars$t1, to=100, by=1)
y <- exp(-pars$r1 * (x - pars$t1))
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")

Close, but not perfect. Here is an example with much smaller time step (0.1) and larger number of cells (10,000).

pars <- list(
  t1 = -60,
  dt = 0,
  r1 = 0.04,
  r2 = 0
)
tchr <- ChromCom3(pars, time = seq(from=-140, to=90, by=0.1))
tchr <- generateCells(tchr, nsim=10000)
x <- seq(from=pars$t1, to=100, by=1)
y <- exp(-pars$r1 * (x - pars$t1))
g <- plotTimelines(tchr)
g + geom_line(data=data.frame(x=x, y=y), aes(x,y), colour="yellow")

Parameter grid

parameterGrid <- function(pars, par1, range1, par2, range2) {
  melts <- NULL
  for(p1 in range1) {
    pars[[par1]] <- p1
    for(p2 in range2) {
      pars[[par2]] <- p2
      chr <- ChromCom3(pars)
      chr <- generateCells(chr)
      label1 <- sprintf("%s = %.3g", par1, p1)
      label2 <- sprintf("%s = %.3g", par2, p2)
      m <- meltTimelines(chr, label1=label1, label2=label2)
      melts <- rbind(melts, m)
    }
  }
  timelinePanel(melts)
}
pars <- list(
  t1 = -30,
  dt = 0,
  r1 = 0.05,
  r2 = 0.01
)
range1 <- c(0.01, 0.05, 0.08, 0.12)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "r1", range1, "r2", range2)

I think we need to introduce a delay between pink and red.

pars <- list(
  t1 = -30,
  dt = 10,
  r1 = 0.05,
  r2 = 0.01
)
range1 <- c(0.01, 0.03, 0.05, 0.10)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "r1", range1, "r2", range2)

What about range of delays?

pars <- list(
  t1 = -30,
  dt = 0,
  r1 = 0.03,
  r2 = 0.01
)
range1 <- c(0, 10, 20, 30)
range2 <- c(0.01, 0.05, 0.08, 0.12)
parameterGrid(pars, "dt", range1, "r2", range2)

Experimental data.

echr <- experimentalData(dataFile$scramble)
plotTimelines(echr, smooth=TRUE, k=15)

13/07/2017

Prepared shiny app for deployment.

24/07/2017

Tomo suggested that the model parameters can be derived directly from data. We can write the following equations

\(\dot{B} = -k_1 B\)

\(\dot{P} = k_1 B - k_2 R\)

\(\dot{R} = k_2 P\)

\(B + P + R = 1\)

However, I’m not sure if this is exacly the model I use in the simulation. In the simulation, we assumed that P->R transition can happen only after the BB->P transition.

Because proportions and their derivatives are observed, we can find the rates:

\(k_1 = - \frac{\dot{B}}{B}\)

\(k_2 = \frac{\dot{R}}{P}\)

This will require a lot of smoothing, otherwise the derivatives will be all over the place.

First, I see if I can recover rates from the generated data.

timeDeriv <- function(chr, k=20) {
  ts <- list()
  for(col in chr$colours) {
    x <- chr$cnt[[col]] / chr$cnt$total
    x <- caTools::runmean(x, k)
    x <- ts(x, start=chr$timepars$start, deltat=chr$timepars$step)
    xdot <- diff(x)
    ts[[col]] <- x
    ts[[paste0(col, ".diff")]] <- xdot
  }
  ts
}
plotRates <- function(chr, k=20) {
  ts <- timeDeriv(chr, k=k)
  k1 <- -ts$BB.diff / ts$BB
  k2 <- ts$R.diff / ts$P
  k1[which(is.nan(k1))] <- NA
  k2[which(is.nan(k2))] <- NA
  df <- data.frame(t=time(k1), k1=k1, k2=k2)
  m <- melt(df, measure.vars = c("k1", "k2"))
  m$value <- as.numeric(m$value)
  m$value[which(m$value > 0.5)] <- 0.5
  m$value[which(m$value < -0.5)] <- -0.5
  ggplot(m, aes(x=t, y=value)) + geom_line(aes(colour=variable)) 
}

Start with a model wiht 10,000 cells and 0.1 min time step. In this model \(k_1 = 0.04\) and \(k_2 = 0\). Top panels show original data, lower panels - smoothed data.

grid.arrange(plotTimelines(tchr), plotRates(tchr, k=1), ncol=2)

grid.arrange(plotTimelines(tchr, smooth=TRUE, k=200), plotRates(tchr, k=200), ncol=2)

What happens if I use only 100 cells and 1 min time step in the model? This time \(k_2 = 0.04\).

chr <- generateCells(chr, nsim=100)
grid.arrange(plotTimelines(chr, k=1), plotRates(chr, k=1), ncol=2)

grid.arrange(plotTimelines(chr, smooth=TRUE, k=20), plotRates(chr, k=20), ncol=2)

And now, for real data. Here are scramble data smoothed with running mean with window size of 10 time points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=10), plotRates(echr, k=10), ncol=2)

Now, smoothing with 80 points.

grid.arrange(plotTimelines(echr, smooth=TRUE, k=80), plotRates(echr, k=80), ncol=2)

LS0tCnRpdGxlOiAiQ2hyb21vc29tZSBjb21wYWN0aW9uIC0gc2ltcGxlIHRocmVlLXN0YXRlIG1vZGVsIgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgICB0b2M6IGZhbHNlCiAgICB0b2NfZmxvYXQ6IGZhbHNlCiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKLS0tCgpgYGB7cn0KbGlicmFyeShDaHJvbUNvbSkKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeShyZXNoYXBlMikKYGBgCgojIyAxMS8wNy8yMDE3CgpXcml0aW5nIGludHJvZHVjdGlvbi4gRmlyc3Qgc2NyaXB0IHRvIGRvIHNpbXVsYXRpb24uCgoKIyMgMTIvMDcvMjAxNwoKTGV0J3MgdHJ5IGlmIGl0IHdvcmtzLiBBbiBleGFtcGxlIGZvciAkdF8xID0gLTYwJCBtaW4sICRcRGVsdGEgdCA9IDAkLCAka18xID0gMC4wNFwge1xybSBtaW59XnstMX0kIGFuZCAka18yID0gMC4wNFwge1xybSBtaW59XnstMX0kLgoKYGBge3IsIGZpZy53aWR0aD01LCBmaWcuaGVpZ2h0PTN9CnBhcnMgPC0gbGlzdCgKICB0MSA9IC02MCwKICBkdCA9IDAsCiAgcjEgPSAwLjA0LAogIHIyID0gMC4wNAopCmNociA8LSBDaHJvbUNvbTMocGFycykKY2hyIDwtIGdlbmVyYXRlQ2VsbHMoY2hyLCBuc2ltPTEwMDApCnBsb3RUaW1lbGluZXMoY2hyKQpgYGAKCkp1c3QgdG8gY2hlY2sgaWYgb3V0cHV0IGlzIGFzIGV4cGVjdGVkLiBXaXRoICRrXzI9MCQgd2Ugc2hvdWxkIGhhdmUgcHVyZSBleHBvbmVudGlhbCBkZWNheS4gVGhlIHllbGxvdyBjdXJ2ZSBpcyAkZV57LWtfMSB0fSQuCgpgYGB7ciwgZmlnLndpZHRoPTV9CnBhcnMgPC0gbGlzdCgKICB0MSA9IC02MCwKICBkdCA9IDAsCiAgcjEgPSAwLjA0LAogIHIyID0gMAopCnRjaHIgPC0gQ2hyb21Db20zKHBhcnMpCnRjaHIgPC0gZ2VuZXJhdGVDZWxscyh0Y2hyLCBuc2ltPTEwMDApCgp4IDwtIHNlcShmcm9tPXBhcnMkdDEsIHRvPTEwMCwgYnk9MSkKeSA8LSBleHAoLXBhcnMkcjEgKiAoeCAtIHBhcnMkdDEpKQpnIDwtIHBsb3RUaW1lbGluZXModGNocikKZyArIGdlb21fbGluZShkYXRhPWRhdGEuZnJhbWUoeD14LCB5PXkpLCBhZXMoeCx5KSwgY29sb3VyPSJ5ZWxsb3ciKQpgYGAKCkNsb3NlLCBidXQgbm90IHBlcmZlY3QuIEhlcmUgaXMgYW4gZXhhbXBsZSB3aXRoIG11Y2ggc21hbGxlciB0aW1lIHN0ZXAgKDAuMSkgYW5kIGxhcmdlciBudW1iZXIgb2YgY2VsbHMgKDEwLDAwMCkuCgpgYGB7ciwgZmlnLndpZHRoPTV9CnBhcnMgPC0gbGlzdCgKICB0MSA9IC02MCwKICBkdCA9IDAsCiAgcjEgPSAwLjA0LAogIHIyID0gMAopCnRjaHIgPC0gQ2hyb21Db20zKHBhcnMsIHRpbWUgPSBzZXEoZnJvbT0tMTQwLCB0bz05MCwgYnk9MC4xKSkKdGNociA8LSBnZW5lcmF0ZUNlbGxzKHRjaHIsIG5zaW09MTAwMDApCgp4IDwtIHNlcShmcm9tPXBhcnMkdDEsIHRvPTEwMCwgYnk9MSkKeSA8LSBleHAoLXBhcnMkcjEgKiAoeCAtIHBhcnMkdDEpKQpnIDwtIHBsb3RUaW1lbGluZXModGNocikKZyArIGdlb21fbGluZShkYXRhPWRhdGEuZnJhbWUoeD14LCB5PXkpLCBhZXMoeCx5KSwgY29sb3VyPSJ5ZWxsb3ciKQpgYGAKCgoKIyMjIFBhcmFtZXRlciBncmlkCgpgYGB7cn0KcGFyYW1ldGVyR3JpZCA8LSBmdW5jdGlvbihwYXJzLCBwYXIxLCByYW5nZTEsIHBhcjIsIHJhbmdlMikgewogIG1lbHRzIDwtIE5VTEwKICBmb3IocDEgaW4gcmFuZ2UxKSB7CiAgICBwYXJzW1twYXIxXV0gPC0gcDEKICAgIGZvcihwMiBpbiByYW5nZTIpIHsKICAgICAgcGFyc1tbcGFyMl1dIDwtIHAyCiAgICAgIGNociA8LSBDaHJvbUNvbTMocGFycykKICAgICAgY2hyIDwtIGdlbmVyYXRlQ2VsbHMoY2hyKQogICAgICBsYWJlbDEgPC0gc3ByaW50ZigiJXMgPSAlLjNnIiwgcGFyMSwgcDEpCiAgICAgIGxhYmVsMiA8LSBzcHJpbnRmKCIlcyA9ICUuM2ciLCBwYXIyLCBwMikKICAgICAgbSA8LSBtZWx0VGltZWxpbmVzKGNociwgbGFiZWwxPWxhYmVsMSwgbGFiZWwyPWxhYmVsMikKICAgICAgbWVsdHMgPC0gcmJpbmQobWVsdHMsIG0pCiAgICB9CiAgfQogIHRpbWVsaW5lUGFuZWwobWVsdHMpCn0KYGBgCgpgYGB7cn0KcGFycyA8LSBsaXN0KAogIHQxID0gLTMwLAogIGR0ID0gMCwKICByMSA9IDAuMDUsCiAgcjIgPSAwLjAxCikKCnJhbmdlMSA8LSBjKDAuMDEsIDAuMDUsIDAuMDgsIDAuMTIpCnJhbmdlMiA8LSBjKDAuMDEsIDAuMDUsIDAuMDgsIDAuMTIpCnBhcmFtZXRlckdyaWQocGFycywgInIxIiwgcmFuZ2UxLCAicjIiLCByYW5nZTIpCmBgYAogSSB0aGluayB3ZSBuZWVkIHRvIGludHJvZHVjZSBhIGRlbGF5IGJldHdlZW4gcGluayBhbmQgcmVkLgogCiAKYGBge3J9CnBhcnMgPC0gbGlzdCgKICB0MSA9IC0zMCwKICBkdCA9IDEwLAogIHIxID0gMC4wNSwKICByMiA9IDAuMDEKKQoKcmFuZ2UxIDwtIGMoMC4wMSwgMC4wMywgMC4wNSwgMC4xMCkKcmFuZ2UyIDwtIGMoMC4wMSwgMC4wNSwgMC4wOCwgMC4xMikKcGFyYW1ldGVyR3JpZChwYXJzLCAicjEiLCByYW5nZTEsICJyMiIsIHJhbmdlMikKYGBgCgpXaGF0IGFib3V0IHJhbmdlIG9mIGRlbGF5cz8KCmBgYHtyfQpwYXJzIDwtIGxpc3QoCiAgdDEgPSAtMzAsCiAgZHQgPSAwLAogIHIxID0gMC4wMywKICByMiA9IDAuMDEKKQoKcmFuZ2UxIDwtIGMoMCwgMTAsIDIwLCAzMCkKcmFuZ2UyIDwtIGMoMC4wMSwgMC4wNSwgMC4wOCwgMC4xMikKcGFyYW1ldGVyR3JpZChwYXJzLCAiZHQiLCByYW5nZTEsICJyMiIsIHJhbmdlMikKYGBgCgpFeHBlcmltZW50YWwgZGF0YS4KCmBgYHtyLCB3YXJuaW5nPUZBTFNFLCBmaWcud2lkdGg9NX0KZWNociA8LSBleHBlcmltZW50YWxEYXRhKGRhdGFGaWxlJHNjcmFtYmxlKQpwbG90VGltZWxpbmVzKGVjaHIsIHNtb290aD1UUlVFLCBrPTE1KQpgYGAKCgojIyAxMy8wNy8yMDE3CgpQcmVwYXJlZCBbc2hpbnkgYXBwXShodHRwczovL3NoaW55LmNvbXBiaW8uZHVuZGVlLmFjLnVrL21hcmVrX2Nocm9tY29tL3BhcmFtX3R1bmVyLykgZm9yIGRlcGxveW1lbnQuCgojIyAyNC8wNy8yMDE3CgpUb21vIHN1Z2dlc3RlZCB0aGF0IHRoZSBtb2RlbCBwYXJhbWV0ZXJzIGNhbiBiZSBkZXJpdmVkIGRpcmVjdGx5IGZyb20gZGF0YS4gV2UgY2FuIHdyaXRlIHRoZSBmb2xsb3dpbmcgZXF1YXRpb25zCgokXGRvdHtCfSA9IC1rXzEgQiQKCiRcZG90e1B9ID0ga18xIEIgLSBrXzIgUiQKCiRcZG90e1J9ID0ga18yIFAkCgokQiArIFAgKyBSID0gMSQKCkhvd2V2ZXIsIEknbSBub3Qgc3VyZSBpZiB0aGlzIGlzIGV4YWNseSB0aGUgbW9kZWwgSSB1c2UgaW4gdGhlIHNpbXVsYXRpb24uIEluIHRoZSBzaW11bGF0aW9uLCB3ZSBhc3N1bWVkIHRoYXQgUC0+UiB0cmFuc2l0aW9uIGNhbiBoYXBwZW4gb25seSBhZnRlciB0aGUgQkItPlAgdHJhbnNpdGlvbi4gCgpCZWNhdXNlIHByb3BvcnRpb25zIGFuZCB0aGVpciBkZXJpdmF0aXZlcyBhcmUgb2JzZXJ2ZWQsIHdlIGNhbiBmaW5kIHRoZSByYXRlczoKCiRrXzEgPSAtIFxmcmFje1xkb3R7Qn19e0J9JAoKJGtfMiA9IFxmcmFje1xkb3R7Un19e1B9JAoKVGhpcyB3aWxsIHJlcXVpcmUgYSBsb3Qgb2Ygc21vb3RoaW5nLCBvdGhlcndpc2UgdGhlIGRlcml2YXRpdmVzIHdpbGwgYmUgYWxsIG92ZXIgdGhlIHBsYWNlLgoKRmlyc3QsIEkgc2VlIGlmIEkgY2FuIHJlY292ZXIgcmF0ZXMgZnJvbSB0aGUgZ2VuZXJhdGVkIGRhdGEuCgpgYGB7cn0KdGltZURlcml2IDwtIGZ1bmN0aW9uKGNociwgaz0yMCkgewogIHRzIDwtIGxpc3QoKQogIGZvcihjb2wgaW4gY2hyJGNvbG91cnMpIHsKICAgIHggPC0gY2hyJGNudFtbY29sXV0gLyBjaHIkY250JHRvdGFsCiAgICB4IDwtIGNhVG9vbHM6OnJ1bm1lYW4oeCwgaykKICAgIHggPC0gdHMoeCwgc3RhcnQ9Y2hyJHRpbWVwYXJzJHN0YXJ0LCBkZWx0YXQ9Y2hyJHRpbWVwYXJzJHN0ZXApCiAgICB4ZG90IDwtIGRpZmYoeCkKICAgIHRzW1tjb2xdXSA8LSB4CiAgICB0c1tbcGFzdGUwKGNvbCwgIi5kaWZmIildXSA8LSB4ZG90CiAgfQogIHRzCn0KYGBgCgpgYGB7cn0KcGxvdFJhdGVzIDwtIGZ1bmN0aW9uKGNociwgaz0yMCkgewogIHRzIDwtIHRpbWVEZXJpdihjaHIsIGs9aykKICBrMSA8LSAtdHMkQkIuZGlmZiAvIHRzJEJCCiAgazIgPC0gdHMkUi5kaWZmIC8gdHMkUAogIGsxW3doaWNoKGlzLm5hbihrMSkpXSA8LSBOQQogIGsyW3doaWNoKGlzLm5hbihrMikpXSA8LSBOQQogIGRmIDwtIGRhdGEuZnJhbWUodD10aW1lKGsxKSwgazE9azEsIGsyPWsyKQogIG0gPC0gbWVsdChkZiwgbWVhc3VyZS52YXJzID0gYygiazEiLCAiazIiKSkKICBtJHZhbHVlIDwtIGFzLm51bWVyaWMobSR2YWx1ZSkKICBtJHZhbHVlW3doaWNoKG0kdmFsdWUgPiAwLjUpXSA8LSAwLjUKICBtJHZhbHVlW3doaWNoKG0kdmFsdWUgPCAtMC41KV0gPC0gLTAuNQogIGdncGxvdChtLCBhZXMoeD10LCB5PXZhbHVlKSkgKyBnZW9tX2xpbmUoYWVzKGNvbG91cj12YXJpYWJsZSkpIAp9CmBgYAoKU3RhcnQgd2l0aCBhIG1vZGVsIHdpaHQgMTAsMDAwIGNlbGxzIGFuZCAwLjEgbWluIHRpbWUgc3RlcC4gSW4gdGhpcyBtb2RlbCAka18xID0gMC4wNCQgYW5kICRrXzIgPSAwJC4gIFRvcCBwYW5lbHMgc2hvdyBvcmlnaW5hbCBkYXRhLCBsb3dlciBwYW5lbHMgLSBzbW9vdGhlZCBkYXRhLgoKYGBge3IsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTN9CmdyaWQuYXJyYW5nZShwbG90VGltZWxpbmVzKHRjaHIpLCBwbG90UmF0ZXModGNociwgaz0xKSwgbmNvbD0yKQpncmlkLmFycmFuZ2UocGxvdFRpbWVsaW5lcyh0Y2hyLCBzbW9vdGg9VFJVRSwgaz0yMDApLCBwbG90UmF0ZXModGNociwgaz0yMDApLCBuY29sPTIpCmBgYAoKCldoYXQgaGFwcGVucyBpZiBJIHVzZSBvbmx5IDEwMCBjZWxscyBhbmQgMSBtaW4gdGltZSBzdGVwIGluIHRoZSBtb2RlbD8gVGhpcyB0aW1lICRrXzIgPSAwLjA0JC4KCmBgYHtyLCBmaWcud2lkdGg9OSwgZmlnLmhlaWdodD0zfQpjaHIgPC0gZ2VuZXJhdGVDZWxscyhjaHIsIG5zaW09MTAwKQpncmlkLmFycmFuZ2UocGxvdFRpbWVsaW5lcyhjaHIsIGs9MSksIHBsb3RSYXRlcyhjaHIsIGs9MSksIG5jb2w9MikKZ3JpZC5hcnJhbmdlKHBsb3RUaW1lbGluZXMoY2hyLCBzbW9vdGg9VFJVRSwgaz0yMCksIHBsb3RSYXRlcyhjaHIsIGs9MjApLCBuY29sPTIpCmBgYAoKCkFuZCBub3csIGZvciByZWFsIGRhdGEuIEhlcmUgYXJlIHNjcmFtYmxlIGRhdGEgc21vb3RoZWQgd2l0aCBydW5uaW5nIG1lYW4gd2l0aCB3aW5kb3cgc2l6ZSBvZiAxMCB0aW1lIHBvaW50cy4KCmBgYHtyLCBmaWcud2lkdGg9OSwgZmlnLmhlaWdodD0zfQpncmlkLmFycmFuZ2UocGxvdFRpbWVsaW5lcyhlY2hyLCBzbW9vdGg9VFJVRSwgaz0xMCksIHBsb3RSYXRlcyhlY2hyLCBrPTEwKSwgbmNvbD0yKQpgYGAKCk5vdywgc21vb3RoaW5nIHdpdGggODAgcG9pbnRzLgoKYGBge3IsIGZpZy53aWR0aD05LCBmaWcuaGVpZ2h0PTN9CmdyaWQuYXJyYW5nZShwbG90VGltZWxpbmVzKGVjaHIsIHNtb290aD1UUlVFLCBrPTgwKSwgcGxvdFJhdGVzKGVjaHIsIGs9ODApLCBuY29sPTIpCmBgYAo=